Comprehensive analyses of solute carrier family members identify SLC12A2 as a novel therapy target for colorectal cancer

As the largest transporter family impacting on tumor genesis and development, the prognostic value of solute carrier (SLC) members has not been elucidated in colorectal cancer (CRC). We aimed to identify a prognostic signature from the SLC members and comprehensively analyze their roles in CRC. Firstly, we downloaded transcriptome data and clinical information of CRC samples from GEO (GSE39582) and TCGA as training and testing dataset, respectively. We extracted the expression matrix of SLC genes and established a prognostic model by univariate and multivariate Cox regression. Afterwards, the low-risk and high-risk group were identified. Then, the differences of prognosis traits, transcriptome features, clinical characteristics, immune infiltration and drug sensitivity between the two groups were explored. Furthermore, molecular subtyping was also implemented by non-negative matrix factorization (NMF). Finally, we studied the expression of the screened SLC genes in CRC tumor tissues and normal tissues as well as investigated the role of SLC12A2 by loss of function and gain of function. As a result, we developed a prognostic risk model based on the screened 6-SLC genes (SLC39A8, SLC2A3, SLC39A13, SLC35B1, SLC4A3, SLC12A2). Both in the training and testing sets, CRC patients in the high-risk group had the poorer prognosis and were in the more advanced pathological stage. What’s more, the high-risk group were enriched with CRC progression signatures and immune infiltration. Two groups showed different drug sensitivity. On the other hand, two distinct subclasses (C1 and C2) were identified based on the 6 SLC genes. CRC patients in the high-risk group and C1 subtype had a worse prognosis. Furthermore, we found and validated that SLC12A2 was steadily upregulated in CRC. A loss-of-function study showed that knockdown of SLC12A2 expression restrained proliferation and stemness of CRC cells while a gain-of-function study showed the contrary results. Hence, we provided a 6-SLC gene signature for prognosis prediction of CRC patients. At the same time, we identified that SLC12A2 could promote tumor progression in CRC, which may serve as a potential therapeutic target.


Western blotting analysis
Proteins were collected from CRC cells and tissues using RIPA lysis buffer (EpiZyme, China), and protein concentrations were detected using a BCA kit (Beyotime, China).30 μg of protein per group was loaded on a 10% SDS-PAGE gel, transferred to a PVDF membrane (Millipore, USA) and then closed with 5% skimmed milk.The blots were cut prior to hybridization with specific antibodies according to the protein molecular weight during www.nature.com/scientificreports/blotting.Thereby, decreasing the amount of incubation solution used during the antibody incubation step, contributing to reduce the use of the costly antibodies.The membranes were incubated overnight at 4 °C in primary antibodies: SLC12A2 (13,884-1-AP, Proteintech, China), C-myc (5605S, CST, USA), Nanog (14,295-1-AP, Proteintech, China), CD44 (37259S, CST, USA), GAPDH (Servicebio, China).Finally, the cells were incubated with the corresponding secondary antibodies for 2 h and the bands were detected by enhanced chemiluminescence reagent (EpiZyme, China).

Cell growth assay
Cell proliferation was detected by Colony Formation and Cell Counting Kit 8 (CCK-8).For Colony Formation assay, 1500 colorectal cancer cells were seeded into 6-well plates with 2 ml of medium per well and cultured for 14 days.The number of colonies with more than 50 cells was counted after 14 days using 4% paraformaldehyde fixation and 0.1% crystal violet staining.As for CCK-8 assay, colorectal cancer cells were inoculated into 96-well plates at a number of 4000 per well and subsequently transfected with siRNA or plasmids for different groups.At specific time points, CCK-8 solution was added to the well plates and cultured for 2 h before detecting the OD450 by means of a microplate reader (BioTek ELx800, USA).

Statistical analysis
Differences in clinical characteristics between the high-risk and low-risk groups were analyzed using chi-square analysis.Student's t-test was used to compare normally distributed variables between the two groups.Two-tailed p-values < 0.05 were statistically significant.

Ethics approval and consent to participate
This project was approved by the Ethics Committee of Zhongnan Hospital of Wuhan University (protocol #2023015 K).Informed consent was obtained from the participated patients.

Identification and validation of a prognosis-related SLC gene signature in CRC
First of all, a flow chart was shown to introduce the design of this study (Fig. 1).The clinical characteristics of GSE39582 and TCGA COAD dataset had been shown in our previous work, and there was no significant difference in general features between two datasets 20 .After excluding the repeated genes and genes whose expression value was zero in any analyzed sample, the expression matrix of the SLC family genes combined with survival data in the two datasets was acquired for subsequent analysis.To screen out prognosis-related SLC genes, univariate cox proportional hazards model was conducted in GSE39582 and TCGA COAD dataset.6 genes (SLC39A8, SLC2A3, SLC39A13, SLC35B1, SLC4A3, SLC12A2) with p-values < 0.05 in the two datasets were regarded as prognosis-related SLC genes.Of the 6 genes, SLC2A3, SLC39A13 and SLC4A3 were risk factors with hazard ratios (HRs) > 1 and SLC39A8, SLC35B1 and SLC12A2 were protective factors with HRs < 1 (Table S1).We also showed the survival curves of the 6 genes (Figure S1).To construct the prognosis-related SLC gene signature, we performed multivariate Cox regression analysis and calculated the risk score of each sample based on the expression and the regression coefficients of the 6 genes (Table S2).Then, patients were assigned to low-risk or high-risk group according to the median risk score in GSE39582 (training set) and TCGA COAD (testing set) dataset.The prognosis of CRC patients in the low-risk group was better than that in the high-risk group in training set and testing set.The heatmap was used to visualize the expression levels of the 6 genes in the low-and high-risk group (Fig. 2A,B).To further identify the importance of the 6-SLC gene signature on the CRC patients' prognosis, we performed survival, ROC and multivariate Cox regression analyses.In the training set, survival analysis revealed that patients with CRC in the low-risk group had a significantly higher survival probability compared to the patients in high-risk group (p = 0.00083) (Fig. 2C-a).ROC analysis showed that the area under the curve (AUC) was 0.62 for the training set (Fig. 2C-b).In addition, multivariate Cox regression analysis also confirmed the independent prognostic value of this signature (Fig. 2C-c).We further applied this signature into testing set and found the consistent results (Fig. 2D).Clearly, these data indicated the 6-SLC gene signature could contribute to prognosis prediction of CRC patients.

Development of nomogram and calibration curves for CRC patients
To accurately estimate the survival of individual CRC patients, we built a nomogram to assess 1-and 3-year survival probabilities based on staging, age, and risk scores in the training and test sets (Fig. 3A,B).Our results suggest that nomograms can be a useful model for prognostic assessment of CRC patients.The calibration curves showed that the predicted prognosis was generally consistent with the actual mortality rates at 1 and 3 years in both the training and test sets (Fig. 3A,B).These data suggest that the nomograms can accurately assess the survival of CRC patients.

Clinical characteristics of CRC patients with different risk
To better clarify the differences of CRC patients in low-risk and high-risk groups, the relationship with several clinical characteristics was studied by Chi-square test.The results in training set was shown in Table 1, which demonstrated that the proportion of patients in "TNM stage", "T stage", "N stage" and "M stage" were significantly different within different risk groups.Consistently, the differences of "TNM stage", "T stage" and "N stage" within different risk groups of testing set also had significance.However, the p-value of "M stage" in testing set was exactly equal to 0.05 (Table 2).A heatmap was further used to more clearly visualize the correlation between different risk groups and clinical characteristics in training set and testing set (Fig. 4A,B).We also explored the correlation between risk groups and "TP53", "KRAS", "BRAF" mutant status in training set.The results showed no significant differences (Table 1).As for MSI status, it showed significant differences in testing set, which was inconsistent with that in training set (Tables 1, 2).Since clinicopathological stage was closely related to the prognosis of patients with CRC, the above results demonstrated the superior performance of the 6-SLC gene signature for prognosis prediction.

Functional enrichment analyses of DEGs in different risk groups
In order to gain insight into the molecular characteristics of low-and high-risk group, we screened the DEGs and conducted their functional enrichment analyses in training dataset.Under a threshold of Adjusted P value < 0.05 and |log2FC|> 0.5, a total of 707 DEGs were identified for the different risk groups.In detail, 536 DEGs were down-regulated while 171 genes were up-regulated (Table S3).GO analysis including biological process (BP), cellular component (CC) and molecular function (MF) showed that the above DEGs were mostly enriched in extracellular matrix (ECM) related functions (Fig. 4C, Figure S2).Moreover, KEGG analysis also exhibited that the pathways of focal adhesion, ECM-receptor interaction, PI3K-Akt and cell cycle were particularly prominent in these DEGs (Fig. 4D).

Immune infiltration and cancer progression estimation in different risk groups
To evaluate the heterogeneity of tumors between low-and high-risk groups, we utilized the ESTIMATE algorithm to calculate stromal scores, immune scores, and tumor purity for both training and test sets.Results showed that the high-risk group had the higher immune score and stromal score, while the tumor purity of the high-risk group was lower than that of the low-risk group (Fig. 5A,B).The training set and the testing set exhibited the highly consistent results.We further investigated the expression of several immune checkpoints in the two risk groups.In training set, the expression of most immune checkpoints in the high-risk group were significantly higher than that in the low-risk group, except for CTLA4 and IL1A (Fig. 5C).Since CD274, CTLA4, IL1A and IL6 were not shown in the filtered expression matrix of the testing set, we analyzed the rest immune checkpoints, the results of which were coincided with that in training set (Fig. 5D).With the significant differences in immune score and immune checkpoints, immune infiltration was investigated to characterize their immunologic landscape.Based on a signature of 17 immune cell type (Table S4), immune cell infiltration was analyzed by GSEA in training set.The heatmap showed that immune cells were significantly enriched in the high-risk group, which was consistent with the result that the group had higher immune scores (Fig. 5E).S5).Results exhibited that the high-risk group displayed higher expression for PI3K-AKT, WNT, MAPK, RAS, NOTCH, HIF-1 and ECM pathways, while the low-risk group was especially enriched with P53 and CELL CYCLE pathways.APOTOSIS pathway had no significance between two groups (Fig. 5F).

Therapeutic drug sensitivity prediction of CRC in different risk groups
Due to the high heterogeneity of CRC, it's important to choose a therapeutic drug with high sensitivity.Therefore, we explored the relationship between the risk score and clinical response to different therapeutic drugs.Based on the Cancer Genome Project (CGP) database, we screened two commonly used chemotherapy drugs (5-Fluorouracil and Cisplatin) and one targeted drug (Cetuximab) for CRC treatment to calculate IC50.The results showed that, both in training set and testing set, CRC patients in the low-risk group were more sensitive to the treatment of 5-Fluorouracil, while patients in the high-risk group had a favorable response to the treatment of Cisplatin.However, there was no significant differences for the treatment of Cetuximab (Fig. 6A,B).Besides, the cancer immunome atlas (TCIA) is a database that provides comprehensive immunogenomic analyses based on the TCGA.Here, we used the TCIA database to evaluate the immunotherapy response of CRC patients with different risk scores through the immunophenoscore (IPS).The results revealed that the total IPS and IPS for CTLA-4 blocker in the low-risk group were significantly higher than that in the high-risk group (Fig. 6C),

Molecular subtyping based on the prognosis-related SLC genes for CRC patients
Recently, molecular subtyping was broadly applied to reveal the tumor heterogeneity.To more comprehensively identify the importance of the 6 prognosis-related SLC genes we selected, NMF clustering was executed.
The consensus map showed that patients with CRC were classified into two distinct clusters in the training set (Fig. 7A-a).There were 344 samples in the cluster 1 (C1) and 212 in the C2.To access the subclasses' assignments, we performed PCA.Data showed that the two clusters were distributed in different side of the two-dimensional coordinate systems (Fig. 7A-b).Survival analysis showed that the survival probability of patients in C1 was significantly lower than that of patients in C2 (Fig. 7A-c).Furthermore, a similar NMF consensus clustering was acquired in the testing set (Fig. 7B-a).Consistently, two subclasses were identified, which also manifested the same distribution as that in training set by PCA (Fig. 7B-b).In the testing set, patients in C1 also had a significantly lower survival probability compared to the patients in C2 (Fig. 7B-c).The correlation analysis of clinical characteristics showed that the proportion of patients in "TNM stage", "T stage" and "N stage" were significantly

SLC12A2 regulated cell proliferation and stemness in CRC cells
The above results demonstrated that SLC12A2 expression was significantly up-regulated in human CRC tissues.
For validation, we collected 29 pairs of tumor tissues and adjacent non-tumor tissues in Zhongnan Hospital of Wuhan University.The results confirmed that SLC12A2 was up-regulated in CRC tumor tissues when compared with non-tumor tissues, both at mRNA level (p = 0.04, n = 29) (Fig. 9A) and protein level (n = 15) (Fig. 9B).Then, SLC12A2 expression level in several kinds of CRC cells was examined.Compared to normal colon cell NCM460, SLC12A2 showed higher expression in CRC cells (Fig. 9C,D).In order to further investigate the biological role of SLC12A2 in CRC cells, two specific SLC12A2 siRNAs (siSLC12A2 #1, #2) were used to knockdown SLC12A2 expression in HT29 and LOVO cells while pCMV3-SLC12A2 plasmids were used to overexpress SLC12A2 expression in HCT116 and SW480 cells.The results of colony formation and CCK-8 assays showed that cell proliferation was distinctly inhibited in SLC12A2 knockdown cells compared with the control cells (Fig. 10A,B).One the other hand, SLC12A2 overexpression induced cell proliferation in HCT116 and SW480 cells (Fig. 10C,D).What's more, western blotting assays showed that knockdown of SLC12A2 down-regulated the expression of the stemness marker c-Myc, Nanog and CD44 while overexpression of SLC12A2 showed the contrary results (Fig. 10E,F).The above data partly indicated that SLC12A2 was involved in CRC progression by promoting cell stemness.

Discussion
As the largest transporter family, SLC membrane transporters had been used to construct the prognosis related signature in lung adenocarcinoma 21 , osteosarcoma 22 , hepatocellular carcinoma 23 , renal cell carcinoma 24 and pancreatic ductal adenocarcinoma 25 .Here, we devoted to reveal a prognosis related SLC gene signature in CRC, which has not been reported before.Our study identified and validated a 6-SLC gene (SLC39A8, SLC2A3, SLC39A13, SLC35B1, SLC4A3, SLC12A2) signature that could predict prognosis of CRC patients.We classified the samples in the training and testing datasets into low-risk or high-risk group based on the expression level of the 6-SLC genes.Different risk group was associated with different prognosis, clinical traits, molecular features, functions, and immune cell fractions.In detail, results showed that CRC patients in the high-risk group had the poorer prognosis.Clinical feature analyses showed that most samples in the high-risk group were in advanced pathological stage.CRC progression signatures, such as PI3K-AKT, WNT, MAPK, RAS, NOTCH, HIF-1 and ECM were also enriched in the high-risk group.Moreover, tumor microenvironment relevant estimation demonstrated that the high-risk group had the higher immune score, stromal score and the lower tumor purity.These data suggested that tumors in the high-risk group were of high heterogeneity and might be refractory.Our opinion was also consistent with the results that the high-risk group had the worse prognosis in both training and testing sets.The high-risk group was filled with a variety of immune cells, and there was a high expression of immune checkpoint genes in the high-risk group, particularly CD274 (also known as PD-L1), CXCR4,CD276, CD4, IL6, LAG3, CCL2, and TGFB1, suggesting that PD-L1 antibodies (e.g., Nivolumab, Durvalumab) and other promising checkpoint inhibitors may be sensitive 26 .Additionally, our analyses revealed that CRC patients in the low-risk group had a www.nature.com/scientificreports/favorable response to the treatment of 5-Fluorouracil and CTLA-4 blocker, which suggested that the 6-SLC gene signature could also be well applied to choose individualized therapies for CRC.
On the other hand, we also developed a molecular subtyping based on the prognosis-related SLC genes for CRC patients.Two distinct subclasses were identified, which showed different prognosis, clinical traits and risk scores.We also showed that patients in the high-risk group and C1 subtype had a worse prognosis.
Furthermore, a detailed analysis of the 6 prognosis-related SLC genes was implemented.We found that SLC35B5 and SLC12A2 were more expressed in malignant cells of CRC compared with other selected prognosisrelated SLC genes.Hence, we evaluated the expression differences of SLC35B5 and SLC12A2 in tumor and normal tissues of CRC.In multiple datasets, we found that SLC12A2 was steadily upregulated in tumor tissues compared with not only over all normal tissues but also paired normal tissues of CRC, which was in consistent with others 27 .What's more, a recent study showed a specific SLC12A2 immunohistochemical staining pattern in precancerous and cancerous colonic UC lesions, which could be helpful for diagnosing dysplasia and cancer in UC and non-UC patients 28 .Thus, we decided to explore the role of SLC12A2 in CRC, which has not been reported before.
The gene SLC12A2 encodes the Na/K/2Cl cotransporter NKCC1, which mainly regulates intracellular ion concentrations and cell volume and plays important functions in neurons and epithelial cells 29 .So far, the role of SLC12A2 in different types of tumors has not been well defined.A study revealed that NKCC1/SLC12A2 was highly expressed in Glioblastomas and it promoted EMT-like process via RhoA and Rac1 signaling pathways 30 .Moreover, inhibition of the NKCC1 could reduce glioma invasion 31 .Detailly, NKCC1 modulated migration of glioma cells by two distinct mechanisms, one was to regulate focal adhesion dynamics and cell contractility; the other was to regulate cell volume through ion transport 32 .Recently, a study showed that blockade of NKCC1  could increase Temozolomide (a conventional chemotherapy drug) induced glioma apoptosis and reduced astrogliosis, which presented the potential therapeutic effect of NKCC1 33 .Besides glioma, it was reported that high expression of NKCC1 predicted poor clinical outcomes for lung adenocarcinoma patients and an EGFR-mutated subgroup 34 .In our study, we proved that SLC12A2 could promote cell growth and stemness in CRC cells.
Here, we presented a pioneer work for identifying a prognosis-related SLC gene signature in CRC.However, we have to mention some defects.Firstly, more datasets are needed to verify our signature.Then, the validation of our signature in clinical samples is necessary.Moreover, our univariate cox proportional hazards model showed that SLC12A2 was a protective factor in CRC while our basic experiments showed that SLC12A2 promoted cancer progression in CRC cells, this inconsistency needs to be explored further.Besides, our data showed that the expression of SLC12A2 was decreased in a stage-dependent manner, the causes also remain unclear.We thought the expression of SLC12A2 might be influenced by the tumor microenvironment, the different therapeutic drugs and so on, which needs to be studied by in vitro and in vivo experiments further.

Conclusions
Overall, our works deepened the understanding of SLC family genes in CRC, and provided a 6-SLC gene signature for prognosis prediction of CRC patients.At the same time, we have tentatively revealed the role of SLC12A2 in CRC, which may serve as a potential therapeutic target.

Figure 1 .
Figure 1.A flow chart of the study.

Figure 2 .
Figure 2. Construction of a prognosis-related SLC gene signature in CRC.(A,B) The distribution of risk scores and the correspondingly survival status of patients, expression abundance of the 6-SLC gene in the training set (A) and the testing set (B). (C) Survival curves(a), ROC curves (b) and multivariate Cox regression analyses (c) in the training set.(D) Survival curves(a), ROC curves (b) and multivariate Cox regression analyses (c) in the testing set.

Figure 3 .
Figure 3. Development of nomogram and calibration curves for CRC patients.(A,B) The nomogram and calibration curves for CRC patient 1-and 3-year survival prediction in the training set (A) and the testing set (B).

Figure 4 .
Figure 4. Clinical characteristics and functional enrichment analyses in different risk groups.(A,B) Heatmap of core clinical characteristics in different risk groups in the training set (A) and the testing set (B). (C) GO analysis for biological process of identified DEGs in the training set.(D) KEGG analysis of identified DEGs in the training set.*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Figure 5 .
Figure 5.Immune infiltration and cancer progression estimation in different risk groups.(A,B) The violin plot of immune score, stromal score and tumor purity from ESTIMATE of two risk groups in the training set (A) and the testing set (B).For violin plots, the three lines within the boxes represent the 25th percentile, median value and the 75th percentile, respectively.The bottom and top of the plots represent the min and max value.(C,D) Expression differences of several immune checkpoint genes between two risk groups in the training set (C) and the testing set (D). (E) Heatmap describing the abundance of immune cell populations in two risk groups.(F) Heatmap describing the abundance of several CRC progression relevant pathways in two risk groups.*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, ns: no significance.

Figure 7 .
Figure 7. Identification of subclasses based on the screened SLC genes for CRC patients.(A,B) Two distinct subclasses (C1, C2) was identified by NMF clustering in the training set (A) and the testing set (B).The consensus map (a), PCA (b) and survival analysis (c) of two subclasses were shown.(C) The risk scores of two subclasses in the training set and testing set.(D) Sankey diagram for the two risk groups and two subtypes.****P < 0.0001.

Figure 8 .
Figure 8. Expression analysis of the 6 prognosis-related SLC genes in CRC.(A) A tSNE map of global cell clusters from 8 CRC samples.(B) Expression distribution in different cell clusters of the 6 SLC genes.(C) Expression levels in malignant cells from 8 CRC samples of the 6 SLC genes.(D) Expression levels of SLC35B5 and SLC12A2 in CRC tumor tissues and normal tissues from GEPIA.(E) Expression levels of SLC35B5 and SLC12A2 in paired tumor tissues and adjacent normal tissues in TCGA COAD and GSE41258.(F) Expression levels of SLC12A2 in the different tumor stage from GEPIA.*P < 0.05.

Figure 9 .Figure 10 .
Figure 9. Validation the expression of SLC12A2 in CRC.(A) The expression of SLC12A2 in 29 pairs of CRC tumor tissues and adjacent non-tumor tissues was tested by real-time PCR.(B) The expression of SLC12A2 in 15 pairs of CRC tumor tissues and adjacent non-tumor tissues was tested by western blotting.(C,D) The expression of SLC12A2 in 6 CRC cell lines and normal colon cells were tested by real-time PCR (C) and western blotting (D).The western blot images in this figure were cropped from full blots and original full blots were presented in Supplementary Fig. 3. Three dependent experiments were performed with similar results.***P < 0.001, ****P < 0.0001.

Table 1 .
Clinical characteristics of patients with different risk in training set.